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Abstract 


The delay Lyapunov equation is an important matrix boundary-value problem 
which arises as an analogue of the Lyapunov equation in the study of time-delay 
systems x(t) = A 0 x(t) + A\x{t — r) + B 0 u(t). We propose a new algorithm for 
the solution of the delay Lyapunov equation. Our method is based on the 
fact that the delay Lyapunov equation can be expressed as a linear system of 
equations, whose unknown is the value U[r/ 2) £ R nxn , i.e., the delay Lyapunov 
matrix at time r/2. This linear matrix equation with n 2 unknowns is solved by 
adapting a preconditioned iterative method such as GMRES. The action of the 
n 2 x n 2 matrix associated to this linear system can be computed by solving a 
coupled matrix initial-value problem. A preconditioner for the iterative method 
is proposed based on solving a T-Sylvester equation MX+X T N = C, for which 
there are methods available in the literature. We prove that the preconditioner is 
effective under certain assumptions. The efficiency of the approach is illustrated 
by applying it to a time-delay system stemming from the discretization of a 
partial differential equation with delay. Approximate solutions to this problem 
can be obtained for problems of size up to n « 1000, i.e., a linear system with 
n 2 « 10 6 unknowns, a dimension which is outside of the capabilities of the other 
existing methods for the delay Lyapunov equation. 

Keywords: Matrix equations, iterative methods, Krylov methods, time-delay 
systems, Sylvester equations, ordinary differential equations 


1. Introduction 

Consider the linear single-delay time-delay system defined by the equations 


(la) 

(lb) 


x(t) = A 0 x(t) + A\x{t — t) + B 0 u(t) 
y{t) = C 0 x(t), 


where Aq,A\ £ R raxn , B 0 £ R nxm , (7)1 £ R nxp . The general equation ([TJ 
appears in many different fields. It is considered a very important topic in the 
field of systems and control, mostly due to the fact that most feedback systems 
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are non-instantaneous in the sense that there is a delay between the observation 
(of for instance the state) and the action of the feedback. See monographs [T] 
and survey paper j3] for literature on time-delay systems. 

The delay Lyapunov equations associated with 0 corresponds to the prob¬ 
lem of finding U £ C°([— t, r],C nxn ) such that 


(2a) 

(2b) 

(2c) 


U\t) 

U(-t) 

—W 


U(t)A 0 + U(t — t > 0, 

U(t ) T , 

U{0)A o + A^U{0) + U(t) t Ai 


A t 


U(t ), 


hold for a given a cost matrix W = W T £ 
instance, W = CqCq). 


(in some applications, for 


Equation (2a) is a matrix delay-differential equation and (2c I is an algebraic 
condition involving C7(0), U{t) and U{—t ) = U(t) t such that ( 2]) can be inter¬ 
preted as a matrix boundary value problem. In this paper we propose a new 
procedure to solve ([2]), with the goal to have good performance for large n. 

The delay Lyapunov equation generalizes the standard Lyapunov equation, 
since e.g., if we set r = 0 the equation reduces to the standard Lyapunov 
equation. Moreover, as established by the last decades of research, the delay 
Lyapunov equation is in many ways playing the same important role for time- 
delay systems as the standard Lyapunov equation plays for standard (delay free) 
linear time-invariant dynamical systems. More precisely, the delay Lyapunov 
equation has been studied in the following ways. It has been extensively used to 
characterize stability of delay differential equations, as one can explicitly con¬ 
struct a Lyapunov functional from 17(f), where the solution is sometimes referred 
to as delay Lyapunov matrices. Sufficient conditions for stability are given in 
a 0 , m and for neutral systems in [7j, and conditions for instability in El mi- 
it has been used to provide bounds on the transient phase of delay-differential 
equations in the PhD thesis Qlj] and [11] [12]. Existence and uniqueness of the 
solutions are well characterized, e.g., in j4j. See also the monograph 2j. Re¬ 
cently, it has been shown that in complete analogy to the standard Lyapunov 
equation the solution to the delay Lyapunov equation explicitly gives the 772- 
norrn m The delay Lyapunov equation can also be used to carry out a model 
order reduction which generalizes balanced truncation m. 

This paper concerns computational aspects of the delay Lyapunov equation. 
Some computational aspects are treated in the literature, e.g., the matrix expo¬ 
nential formula in uni, the polynomial approximation approach in HE spectral 
(Chebyshev-based) discretization approaches in [TSJ 1(1 and an ODE-approach 
in the PhD thesis [T7i Chapter 3]. 

In complete contrast to the delay Lyapunov equation, the computational 
aspects of the standard Lyapunov equation have received considerable atten¬ 
tion, mostly in the numerical linear algebra community. Most importantly, the 
Bartels-Stewart method [18] . ADI methods m , Krylov methods mm, and 
rational Krylov methods [22j . including preconditioning techniques [23]> have 
turned to be effective in various situations. For a more thorough review, see 
the survey [24]. To our knowledge, there exist no natural generalization of the 
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U(~T ) 


symmetry 


U{t/2) U(t) 

• =? 

Z 2 (0) ^(0) Zi(t/2) 


U( 0 ) 


W2) 


Figure 1: Graphical representation of the relation between U(t), Z\(t) and Z 2 (t). 


Bartels-Stewart algorithm and there are no Krylov methods for delay Lyapunov 
equation. 

The method we propose is tailored to medium-scale equations; it combines 
the use of a Krylov-type method and a direct algorithm similar to the Bartels- 
Stewart one. More precisely, our approach is based on a characterization of the 
solution to the delay Lyapunov equation as a linear system of equations with 
n 2 unknowns. This characterization is derived in Section [U Since the linear 
system derived in Section [2] is large and only given implicitly as a matrix vector 
product, we propose to adapt iterative methods which are based on matrix vec¬ 
tor products only, e.g., GMRES [23] or BiCGStab [2H] , to this problem. It turns 
out to be natural to use a preconditioner involving a matrix equation called the 
T-Sylvester equation, for which there are efficient 0(n 3 ) methods for the dense 
case EZj. We quantify the quality of the preconditioner by deriving a bound 
on the convergence factor of the iterative method. The iterative method and 
the preconditioner are given in Section [3] The performance of the approach is 
illustrated with simulations in Section ffl We apply the method to a problem 
stemming from the discretization of a two-dimensional partial delay-differential 
equation (PDDE). The number of iterations appears to be essentially indepen¬ 
dent of the grid, which suggests that the preconditioner is a sensible choice for 
this PDDE. 

We use notation which is standard for analysis of matrix equations. The 
vectorization operation is denoted vec(f?), i.e., if B = [&i ... & m ] G R nxm , 

vec {B) T = [frf ... . The Kronecker product is denoted (g). LTnless other¬ 

wise stated, || ■ || denotes the Euclidean norm for vectors and the spectral norm 
for matrices. We denote the Frobenius norm by || • ||^. 

2. Reformulation of the delay Lyapunov equations 

Our method is based on a reformulation of the delay Lyapunov equation 
where we define for each t G [0, r/2] 

(3) Z x (t) :=U(T/2 + t), Z 2 (t):=U{T/2-t). 

The two matrix-valued functions Zi(t) and Z 2 (t) coincide with U(t) up to a 
change of the time coordinate which is represented visually in Figure [l] Es¬ 
sentially, they represent two different branches of U(t) “taking off” from t/ 2 
in opposite directions. Note that the left half of the function, U([— r,0]), is 
determined uniquely by the right half f/([0,r]) by the transposition symmetry 
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condition (2b I. The only nontrivial condition implied by (2b) is that U (0) must 
be symmetric. 

Note that 

(4&) Zi(t — t) = U(t — t + t/ 2) = U(t — t/ 2) = U(r/2 — t) T = Z 2 (t) T 
(4b )Z 2 {t-r) = Ufir/2-t-r) = U(-t-r/2) = U{t + T/2) T = Zf(t) 


Hence, the delay differential equation (2a) becomes an ordinary differential equa¬ 
tion 


(5a) Z' x (t) — Zi(t)Ao + Z 2 (t) T Ai, 

(5b) Z' 2 (t) = -Z 1 (t) T A 1 - Z 2 (t)A 0 . 


This is a constant-coefficient homogeneous linear system of ODEs which can be 
solved explicitly if the common (unknown) initial value Z x (0) = Z 2 { 0) = U(r/2) 
is provided. Using vectorization, we can give an explicit formula 


( 6 ) 


vec Zi(t) 
vec Z 2 \t) T 


exp(L4) 


vec U(t/ 2) 
vec U(t/2) t 


where 


( 7 ) 


Aq In Af (§) In 

—In ® Af —I n ® A^ 


In terms of Z\ (t) and Z 2 (t ), the algebraic condition (2c) and the symmetry 
condition (2b) for t = 0 reduce to 


(8a) 0 = W + Z 2 {t/2) t A 0 + A^Z 2 {t/2) + Z 1 {t/2) t A 1 +AIZ 1 {t/2 ), 
(8b) 0 = Z 2 (t/2) — Z 2 [t/ 2) t . 


Notice that the right-hand side of (8a) is symmetric and that of (8b) is anti¬ 
symmetric. A linear combination of them gives 


(9) 0 = W+Z 2 (t/2) t (A 0 - cl) + (Aq + cI)Z 2 (t/2) + Z x (r/2) T A x +Af Z i (r/2) 


for each c £ R, which forms the basis of our matrix operator. 
Definition 1 . Let L c : R nxrl —> R raxn be defined by 


(10) L C (X)~ 

Z 2 (t/2) t (Aq - cl) + « + cI)Z 2 (t/ 2 ) + Z 1 {t/2) t A 1 + A x Z x {t/ 2 ) 

where Zi : [0, r/2] —» K nxn , i = 1,2 are the unique solutions to the initial value 
problem ([5]) with Z x (0) = Z 2 { 0) = X. 

We shall need the following easy linear algebra result. 

Lemma 2. Let M = M T £ R nxrl and N = —N T £ M nxn be two matrices, one 
symmetric and one antisymmetric. Then, M+N = 0 if and only if M = N = 0. 
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Proof. The ‘if’ part is trivial; let us prove the ‘only if’. Suppose M + N = 0; 
then, transposing, we have also 0 = M T + N T = M — N. Summing and 
subtracting the two relations we have 2 M = 2TV = 0. □ 

A time-delay system is called exponentially stable if ||x(i)|| < Mexp(-fit) 
for some constants M > 0, /3 > 0. If this condition holds, then the solution U ( t ) 
to 0 is unique da Theorem 4], In this case, we can formulate the equivalence 
between the delay Lyapunov equation and a linear system with operator L c . 

Theorem 3 (Equivalence). Suppose Aq and A\ and r are such that |l]) is 
exponentially stable and let W £ R nxn be any symmetric matrix. Let U be the 
solution to the delay Lyapunov matrices ([2| and let L c be defined by ( |10| ) . Then, 
for any c 0, X = U(t/ 2) is the unique solution of the linear system 


(11) L c { X) = — W. 

Proof. Equation ([9]) already shows that if X = U(t/ 2) then L C (X) = —W. It 
remains to prove the reverse implication. Suppose that X satisfies L C (X) + W = 
0; then, by Lemma [2] applied to 

M = Z 2 (t/2) t A 0 + Aq Z 2 (t /2) + Z 1 (t/2) t A 1 + A^Zfir/ 2) - W, 

N = c(Z 2 (t/2) - Z 2 (t/ 2) t ), 


the conditions <[8j) hold. Define 


Z 2 (r/2 — t) 0 < t < t/2, 
U(t)= ^ Z\(t — t/2) t/2 < t < r, 
U{—t) T ~T<t< 0. 


The function U(t) is continuous in 0 by (8b), and in ±r/2 by the choice of initial 
conditions, hence it is globally continuous on [—r, t). Moreover, the differential 
equation (2a) holds for all t ^ 0,r/2. By continuity, it must also hold for these 
values. Hence U(t) solves ([ 2 ]). As we assume exponential stability, the solution 
is unique and hence U(t) = U(t). □ 


Since the linear system L C (X) = — W has a unique solution for each sym¬ 
metric W £ R" xn , we have the following result. 

Corollary 4. Suppose 0 is exponentially stable. Then, the linear operator L c 
is nonsingular for each c / 0. 


A delay-free formulation of the delay Lyapunov equations has also been 
derived in [U Equation (13)]. That formulation cannot be described with a 
linear operator in a way that can be adapted to an iterative method in the same 
way that we show in the following section. 
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3. Algorithm 


We now know from the previous section that the matrix equation © is 
equivalent to the delay Lyapunov equation. By vectorizing (11), we obtain the 
linear system on standard form 


( 12 ) 


vec L c (ve c 1 x) = — vec W, 


where the inverse function vec : (a;) maps vec A £ K" to A £ 


Let 


A £ R" xn the matrix associated to it. We know that A is nonsingular by 
Corollary [4] 

Our approach is based on specializing an iterative method for linear systems 
to ©■ In order to specialize an iterative method for large-scale linear systems, 
we need two ingredients. We need an efficient procedure to compute the action 


corresponding to the left-hand side of (12); and we need a preconditioner. These 


two ingredients are described in the following two subsections. 


3.1. Action of L c 


The action of the operator L c is defined by ([5]) and (101. As a consequence, 
the recipe to compute L C (A) for a given matrix X is simple: 


1 . 

2 . 


Compute the solutions Z\{t/2), Z 2 {t/2) of the linear, constant-coefficient 
initial-value problem © with initial values Zi(0) = Z 2 (fi) = X. 

Compute L C (X) using the expression (10). 


In practice, a detail is crucial in the choice of the numerical algorithm for the 
first step. We distinguish two possible scenarios: 


• We use a method with a fixed step-size and no adaptivity: for instance, 
the (explicit or implicit) Euler method, or a non-adaptive Runge-Kutta 
method. In this case, we are effectively substituting L c with a different 
operator L c , which replaces the differential operator in Step [l] with a finite 
discretization. This operator (for most classical methods) is still linear, so 
the theory of Krylov subspace methods can be applied without changes: 
we are applying a Krylov method to get an approximate solution of a 
nearby linear problem L c . 


• We use an adaptive method, which can change step size along the algo¬ 
rithm, possibly in different ways for different initial values A'. For instance, 
the Runge-Kutta-Fehlberg method (Matlab’s ode45). While apparently 
the two cases are similar, the addition of adaptivity has an important 
consequence: the computed operator L c , this time, is no longer a linear 
operator, because in general L C {X\ + A 2 ) / L c (Xi) + L C (X 2 ). Indeed, for 
different values of the input A the initial-value problems could be solved 
using different grids, and hence different discrete approximations of the 
propagation operator. Thus we are in the realm of inexact Krylov meth¬ 
ods |28j . The theory in this case is more involved, as more care is required 
with the error thresholds. 
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It turns out in the simulations (in Section [4]) that it is advantageous to try to 
keep the number of iterations as low as possible. In this paper we therefore 
focus on the first approach, since an inexact Krylov method can require more 
iterations. 


3.2. Preconditioning 

In order to make iterative methods effective, it is common to carry out a 
transformation which preconditions the problem. This can often be interpreted 
as transforming the problem with an approximation of the inverse of the ma¬ 
trix/operator. We focus on a particular preconditioner obtained by solving the 


problem exactly when A\ is replaced with the zero matrix. Then (10) becomes 


(13) 


L C (X) := Z 2 (t/2) t (A 0 - cl) + (A% + cI)Z 2 {t/2 ), 


and (5b) decouples from Z\ such that Z' 2 = —Z 2 (t)Ao, which we can solve 
explicitly to get Z 2 (r/ 2) = X exp(— tAq/2). 

Let T be the operator 


T(Y) = (A% + cI)Y + T t (H 0 - cl). 


The operator L c is invertible if and only T 1 exists, and in this case we have 
(14) L-\Z)=T-\Z)eMrA 0 /2). 


Inverting the operator T correspond to solving the so-called (real) T-Sylvester 
equation MY + Y T N = C. The paper m discusses the solvability of this equa¬ 
tion and presents a direct 0(n 3 ) Bartels-Stewart-like algorithm for its solution. 
In particular, the following result holds. 

Theorem 5 ([29l Lemma 8 ], [27] ). Let M,N,C € R nx ”. The equation MX + 
X T N = C has a unique solution X for each right-hand side C if and only if 
fiijlj 7 ^ 1 for each pair pi, p,j of eigenvalues of the pencil M — XN T . 

In our case, M = A^ + cl, N = A 0 — cl, so after a quick computation the 
solvability condition reduces to the following condition, which is independent of 

c. 


Definition 6 (Hamiltonian eigenpairing). We say that the matrix Ho £ K raxn 
has no Hamiltonian eigenpairing, if for each pair of eigenvalues \i,\j of the 
matrix Aq, we have 

Xi + A j 0. 

A matrix has no Hamiltonian eigenpairing, for instance, if < 0 for each 
eigenvalue A of H 0 , i.e., if the delay-free system obtained by setting Hi = 0 is 
stable. 

In order to characterize the convergence and quality of the preconditioner 
we use a fundamental min-max bound. Suppose we carry out GMRES on the 
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matrix A £ M. NxN with eigenvalues Ai,...,Ajv- From [25], Proposition 4] we 
have the bound of the residual 

\\r m+1 \\< K (X)e^\\r 0 \\ 


and 


g( m ) = min max|p(Aj)| 

p£Pm i 


where P m = {p : polynomial of degree m such that p(0) = 1}. We now apply 
the standard Zarantonello bound EM Lemma 6.26], where we assume that the 
eigenvalues are contained in a disk of radius r centered at c = 1, correspond¬ 
ing to selecting p(z) = such that < r m /c m = r m < ||A — I\\ m . 

Since preconditioned GMRES with preconditioner A is essentially equivalent 
to GMRES applied to the matrix A~ X A, a bound on ||T _1 Al — /|| provides a 
characterization of the convergence factor of preconditioned GMRES. Because 
of the vectorization included in our setting, bounding ||Al _1 A — I1 | corresponds 
to giving an estimate for the quantity 


\\L~\L C (X)) - X\\ F 
\\X\\f 


Our preconditioner is constructed by setting A\ = 0. Therefore, we expect that 
the preconditioner works well if ||Ai|| is small. This reasoning is formalized in 
the following result. 


Theorem 7 (Quality of preconditioner). Suppose the system Q is exponen¬ 
tially stable and suppose that Aq has no Hamiltonian eigenpairing. Let L c and 
L c be defined by (10) and © respectively. Then, 


(15) 


\\L- C \L C (X)) - X\\ F 

imif 


o(Piii 2 ), 


where the constant hidden in the O(-) notation depends only on ||Ro||, r and c. 


Proof. We invoke Lemma [9] to bound the left-hand side of (15) 


(16) 


L- 1 (L C (X)) - X 


" -1 


(l c (X) - L c (xj) 


I All 


|A|| 


K exp(r||A 0 ||/2) 


< 


L C (X)-L C (X) 


\X\\ 


In order to bound L C (X) — L C {X) we let Z\ and Z 2 correspond to L C (X), 
i.e., they satisfy the equations ([5]) with initial value Z i(0) = Z 2 (0) = X. Let Z\ 
and Z 2 be the corresponding quantities for L C {X). Moreover, let A 2 := Z 2 — Z 2 . 
We have 


(17) L C {X)-L C {X) = 

A 2 (t/2) t (A 0 - cl) + {A? + c/)A 2 (r/2) + Z 1 (t/2) t A 1 + A^Z^t/2), 























for which A 2 (t/ 2) and Z\ (t/ 2) can be bounded as follows. Lemma [ 8 ] tells us 
that 

(18) \\Z 1 (t/2)\\ f < 2exp(T(||A 0 ||2 + ||Hi|| 2 ))||X|| f . 

By definition, A 2 satisfies the ODE 

(19) A' 2 {t) = -A 2 {t)A 0 + g{t), A 2 (0) = 0, 


where g(t) := — Zi(t) T A\. The variation-of-constants formula applied to (19) 
results in the explicit expression 

r t 


A 2 (i) = - / Zi{s) T Ai exp((s - t)A 0 ) ds. 
Jo 


Hence, 


,7-/2 

(20a) ||A 2 (t/2)|| f < / \\Zi(s) T A\ exp((s — t/2)Aq)\\ f ds 

Jo 

/2 


( 20 b) 

( 20 c) 


< 


|Zi(s)llFPill 2 ll ex P(( s - r/2)A 0 )\\ 2 ds 


/o 


< rexp(r(||H 0 || 2 + ||A 1 || 2 ))||H 1 || 2 exp(r||A 0 || 2 / 2 )||X|| F . 


We now evaluate the Frobenius norm of (17) and apply the triangle inequality 
and the bounds (18) and ( 20 ), which shows that 

( 21 ) 


\L C (X) - L c (X)\\ f 


II All 


= 0(Pi|| 2 ). 


The hidden constant in (21) depends only on ||Ao|| 2 , c, and r. The conclusion 
(15) follows by combining (16) and (21). □ 


4. Simulations 


4-1. A small example 

In order to illustrate the preconditioner and properties of our approach we 
first consider a small example with randomly generated Aq matrix. We specify 
the matrices for reproducibility 



-26 

22 

-1 

-4 

2 

-24 

-4 

1 

7 

11 

-24 

-22 

-13 

15 

-1 

-9 


Hi = odiag(— 1 , —0.5,0,0.5), W = I 


and t = 1. We carry out simulations for different a = ||^4i||. The time-delay 
system is stable for all a € [0,10]. The corresponding delay Lyapunov equation 
satisfies 


U(t/ 2 ) 


1 

loo 


0.2302 -0.0156 

-0.0885 0.0044 

0.1466 -0.0057 

-0.5485 0.0331 


0.0101 -0.3729 

-0.0038 0.1380 

0.0056 -0.2263 

-0.0238 0.8755 


9 




















for a = 1. 

We combine our approach with two different generic iterative methods for 
linear systems of equations, GMRES j25| and BiCGStab [26] and select c = 1. 
To illustrate the properties of the performance of the iterative method, we solve 
the ODE defining L c to full precision with the matrix exponential. The absolute 
error as a function of iteration is given in Figure [2] Both methods successfully 
solve the problem before the break-down at iteration n 2 except for ||Ai|| = 10. 
No substantial difference between the two iterative methods can be observed 
in the error as a function of iteration, i.e., nothing can be concluded regarding 
which of the two variants is better for this problem. The convergence of the two 
methods is faster for small ||Ai||. This is due to the fact that the preconditioner 
is more effective when ||Ai|| is small, which is consistent with Theorem [T] and 
Figure [3] where we clearly see that the norm of the preconditioned system X 1 —»■ 
L~ 1 (L C (X)) has a linear dependence on ||Ai||. The same conclusion is supported 
by the localization of the eigenvalues of the linear map X 1 —>■ L~ 1 (L C (X)) in 
Figure [3] 0 . 


Et, 


'X 

I 

* 


3 

O 

co 




10 2 

10 ° 

10 -2 

10 -4 

10 -6 

10 -8 

10-io 

10- 12 

10 -14 



Iteration k 


Figure 2: Convergence for different preconditioned iterative methods applied to the small 
example in Section 1 4. 1| 


4-2. A large-scale example 

In relation to other methods for delay Lyapunov equations, our iterative ap¬ 
proach is likely to have better relative performance for large problems. We illus¬ 
trate this with the following time-delay system stemming from the discretization 
of a partial differential equation with delajQ More precisely, we consider on the 


1 The MATLAB-code for the example and the simulation is publicly available on http: 

//www.math.kth.se/~eliasj/src/dlyap_precond 
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(a) The preconditioner is better for small ||Ai|| as 
predicted by Theorem [?] 


• 10~ 5 - 10~ 3 10~ 2 


5 
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4 
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1 
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1,000 
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X 


o 

2 

- 

0.5 

— 



0 

* X - 

0 

§0 

0 

- 

0 

- o- 

0 

• • 

5 

X 

T X _E 

-1 

o 

1 Q _L_ 

-2 

-4 

-i q_^ 

-0.5 

-1 

_j_= 

-1,000 

• 


0.9999 1 1.0001 0.995 1 1.005 0.95 1 1.05 0 1 2 3 -1,000 0 1,000 


(b) Real and imaginary parts of the eigenvalues of the operator X > L c ( L C {X )) for different 
a with symbols specified in Figure [2] 


Figure 3: Illustration of the quality of the preconditioner. 


domain (x, y) G [0,1] x [0,1] the PDDE 

dv 

( 22 a) v(x,y,t) = Av(x,y,t) + v(x,y,t) + f{x, z) — (x,y,t - t) + u(t.) 
( 22 b) w(t) = v(l/2, 1 / 2 ) 


where f(x , y) = fo cos (xy) sin( 7 rx) with homogeneous Dirichlet boundary condi¬ 
tions, and fo = 5. The PDDE (22) can be interpreted as waves propagating on 
a square, with damping and delayed feedback control. PDDEs are for instance 
studied in m, In order to reach a problem of the form ([T]) we rephrase ( [22] ) as a 
system of PDDEs which is first-order in time. We carry out a semi-discretization 
with finite differences in space with n x +1 intervals in the x-direction and n y +1 
intervals in the y-direction, i.e., h x = l/{n x + 1 ), x = kh x , k = 1 ,... ,n x and 
h y = 1 /(n y + 1), yk = kh y , k = 1 ,...,n y . The corresponding discretized 
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An — 


I 0 D XX T Dyy 


I 

>7 —7 


time-delay system is of the form (JTJ) with coefficient matrices given by 
(23a) J ° 

(23b) 

(23c) 


Ai — 

B 0 = [1 


0 O' 

diag(7 ;l )(J (g> D x ) 0 

• 1 0 ••• 


of 


(23d) 


where 


D x x 


Dx = 


1 

hi 


1 

2 hx 


Co = 


r -2 l 


l ' • ‘ • 


e (n B +l)/2 ® e (n x +l)/2 0 


1 - 2-1 


% Dyy ~ hi 


r -2 l 


l * • ’ • 


' • * • l 

1 - 2-1 


0 1 


-1 ’• *• 


* • ‘ • 1 
-1 0 


F = vec([/(®i,i/j)]"j=i)- 


In the setting of T^morm computation (as in E9) we need to solve the delay 
Lyapunov equation with W = CqCq. 

We carried out simulations of this system using a computer with an Intel i7 
quad-core processor with 2.1GHz and 16 GB of RAM. For the finest discretiza¬ 
tion that we could treat with our approach, we have n x = n y = 23, n = 1058, 
||^4 0 II 2 ~ 5000 and ||Alj || « 100. We again select c = 1. 

In order to solve the ODE ([5]) we used a fixed fourth order Runge-Kutta 
method with N = 500 grid points. The iteration history of the two variants 
is visualized in Figure [4] for n = 1058. We observe linear convergence and no 
substantial difference in convergence rate. The execution time of our approach, 
in relation to some other approaches in the literature are reported in Table [T] 
We observe better relative computation time for larger problems. Moreover, we 
observe that other approaches fail due to requirement on memory resources. 

Note also in Table |T] that the number of iterations required to reach a 
specified tolerance appears not to grow substantially with the size of problem. 
Hence, the method appears to have essentially grid-independent convergence 
rate, which is considered a very important feature of a preconditioner. 

In a detailed profiling of our approach, we identify that two components are 
dominating, solving the ODE, i.e., computing the action, and solving the T- 
Sylvester equation. For the finest discretization, solving one T-Sylvester equa¬ 
tion took approximately 320 seconds and carrying out one step of the ODE 
required 30 seconds. Since the main computational effort lies on the solution of 
the T-Sylvester equation and not in the computation of the action of L c , it is 


for this problem no advantage to use adaptivity (as discussed in Section 3.1). In 
fact, the number of iterations can even increase when inexact Krylov methods 
are used. 
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We note that the implementation that we have used to solve T-Sylvester 
equations is not particularly optimized; it is a vectorized version of the algorithm 
in m that we have implemented in MATLAB for use in these experiments. The 
complexity in flops of the required computations is only slightly larger than what 
is required for solving a standard Sylvester equation with the Bartels-Stewart 
algorithm, a task which requires less than 8 seconds on our machine. Hence, we 
expect a major reduction in the timings if a carefully optimized solver for the 
T-Sylvester is used instead. 

To our knowledge, the largest delay Lyapunov equation previously solved is 
with n = 110 in [T3] , 



Figure 4: Convergence of the iterative methods with T-Sylvester preconditioning corre¬ 
sponding to the time-delay system stemming from the discretization of the PDDE \22\ with 
n = 2n x n y = 1058 for the example in Section |4.2| 



Matrix exp. [10] 

Discr. first 

Proposed method 


Wall time 

Wall time 

Wall time 

iterations 

n = 28 

0.5 s 

0.06 s 

3.9 s 

19 

n = 50 

296.4 s 

0.6 s 

12.2 s 

21 

n = 242 

MEMERR 

30.9 s 

232.8 s 

25 

n = 722 

MEMERR 

MEMERR 

4629.3 s 

27 

n = 1058 

MEMERR 

MEMERR 

3.1 hours 

27 


Table 1: Performance in relation to other methods. MEMERR represents runs which could 
not solve the problem to reasonable accuracy. Our approach was run with GMRES, RK4 (with 
N = 500) and termination tolerance 10 —12 . Discr. first represents the approach discussed in 
m and used in El with N = 10 grid points, which resulted in a solution of accuracy in U( 0) 
in general larger than 10 -7 . 
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5. Concluding remarks and outlook 


We have in this paper proposed a procedure to solve delay Lyapunov equa¬ 
tions based on iterative methods for linear systems combined with a direct 
method for T-Sylvester equations. Although the method performs well in prac¬ 
tice, there appears to be possibilities to improve it further, which we consider 
beyond the scope of the paper. 

As observed in the simulations, the dominating ingredient of the approach 
is the solution to the T-Sylvester equation. Hence, in order to solve even larger 
problems we need new methods for T-Sylvester equations. Improvements are 
possible, e.g., by lower level implementations, or by developing methods which 
can take the sparsity of the matrices into account, e.g., similar to the Krylov 
methods and rational Krylov methods for Lyapunov equations [2D] or approaches 
based on Riemannian optimization . ’52] , 

The preconditioner in general plays an important role in iterative methods 
for linear systems and the effectiveness of the preconditioner is typically very 
problem-dependent. This is also the case in our approach. Although the sim¬ 
ulations often worked well, during some experiments, in particular situations 
where Ao have some eigenvalues which are very negative, the preconditioner did 
not appear very effective, even if ||Ai|| was quite small. This can be due to the 
fact that the hidden constant in the expression (15) may be large. 

Since the action of L c is based on solving an initial value problem, there 
are very natural options to carry out action of L c in an inexact way, with 
a predefined tolerance, thereby saving computation cost. There exists theory 
for inexact Krylov methods, in e.g., [28] . which could be applied in such an 
approach. However, in our setting, where computing the action of the precon¬ 
ditioner is a dominating cost, an approach with inexact matrix vector product 
would require more iterations, i.e., more preconditioned actions, and does not 
appear advantageous in this setting. 

The delay Lyapunov equation has been generalized in several ways, e.g., to 
multiple delays and neutral systems. Our approach might be generalizable to 
some of these cases. The simplest situations appears to be if the delays are 
integer multiplies of each other, also known as commensurate delays. For the 
commensurate case there are procedures which resemble our reformulation ff 
with Sylvester resultant matrices cm Problem 6.72], However, this increases 
the size of the problem. An attractive feature of our approach is that we work 
only with matrices of size n, which would not be the case in the direct adaption 
to multiple commensurate delays using [101 Problem 6.72]. 


Acknowledgments 

The authors thank Antti Koskela and Tobias Damm for discussions about 
early results of the paper. 

F. Poloni acknowledges the support of the PRA 2014 project “Mathematical 
models and computational methods for complex networks” of the University of 
Pisa, and of INDAM (Istituto Nazionale di Alta Matematica). E. Jarlebring 


14 



acknowledges the support of the Swedish research council (Vetenskapsradet) 

project 2013-4640. 

References 

References 

[1] W. Michiels, S.-I. Niculescu, Stability and Stabilization of Time-Delay Sys¬ 
tems: An Eigenvalue-Based Approach, Advances in Design and Control 12, 
SIAM Publications, Philadelphia, 2007. 

[2] K. Gu, V. Kharitonov, J. Chen, Stability of Time-Delay Systems, Control 
Engineering. Boston, MA: Birkhauser, 2003. 

[3] J.-P. Richard, Time-delay systems: an overview of some recent advances 
and open problems, Automatica 39 (10) (2003) 1667 -1694. 

[4] V. Kharitonov, Lyapunov matrices for a class of time delay systems, Syst. 
Control Lett. 55 (7) (2006) 610-617. 

[5] G. Ochoa, V. Kharitonov, S. Mondie, Critical frequencies and parameters 
for linear delay systems: A Lyapunov matrix approach, Syst. Control Lett. 
26 (2013) 781 -790. 

[6] G. Ochoa, V. Kharitonov, Lyapunov matrices for neutral type time delay 
systems, in: Proceedings of the 2nd International Conference on Electrical 
and Electronics Engineering, Mexico City, Mexico, 2005. 

[7] G. Ochoa, J. Velazquez-Velazquez, V. Kharitonov, S. Mondie, Lyapunov 
matrices for neutral type time delay systems, in: Proceedings of the 7th 
IFAC workshop on time delay systems, Nantes, France, 2007. 

[8] S. Mondie, G. Ochoa, B. Ochoa, Instability conditions for linear time delay 
systems: a Lyapunov matrix function approach, Int. J. Control 84 (10) 
(2011) 1601-1611. 

[9] A. Egorov, S. Mondie, Necessary stability conditions for linear delay sys¬ 
tems, Automatica 50 (12) (2014) 3204-3208. 

[10] E. Plischke, Transient effects of linear dynamical systems, Ph.D. thesis, 
Universitat Bremen (2005). 

[11] V. Kharitonov, D. Hinriclrsen, Exponential estimates for time delay sys¬ 
tems, Syst. Control Lett. 53 (5) (2004) 395-405. 

[12] V. Kharitonov, E. Plischke, Lyapunov matrices for time-delay systems, 
Syst. Control Lett. 55 (9) (2006) 697-706. 

[13] E. Jarlebring, J. Vanbiervliet, W. Michiels, Characterizing and computing 
the T-L -2 norm of time-delay systems by solving the delay Lyapunov equation, 
IEEE Trans. Autom. Control 56 (4) (2011) 814-825. 


15 



[14] E. Jarlebring, T. Damm, W. Michiels, Model reduction of time-delay sys¬ 
tems using position balancing and delay Lyapunov equations, Math. Con¬ 
trol Signals Syst. 25 (2) (2013) 147 -166. 

[15] E. Huesca, S. Mondie, J. Santos, Polynomial approximations of the Lya¬ 
punov matrix of a class of time delay systems, in: Proceedings of the 8th 
IFAC workshop on time-delay systems, Sinaia, Romania, 2009. 

[16] J. Vanbiervliet, W. Michiels, E. Jarlebring, Using spectral discretisation for 
the optimal W? design of time-delay systems, Int. J. Control 84 (2) (2011) 
228-241. 

[17] A. Merz, Computation of generalized gramians for model reduction of bi¬ 
linear control systems and time-delay systems, Ph.D. thesis, T.U. Kaiser¬ 
slautern (2012). 

[18] R. Bartels, G. W. Stewart, Solution of the matrix equation AX + XB = C, 
Comm A.C.M. 15 (9) (1972) 820-826. 

[19] P. Benner, J.-R. Li, T. Penzl, Numerical solution of large-scale Lya¬ 
punov equations, Riccati equations, and linear-quadratic optimal control 
problems, Numer. Linear Algebra Appl. 15 (9) (2008) 755-777. doi: 
10.1002/nla.622 

URL http://dx.doi.org/10.1002/nla.622 

[20] V. Simoncini, A new iterative method for solving large-scale Lyapunov 
matrix equations, SIAM J. Sci. Comput. 29 (3) (2007) 1268-1288. 

[21] D. Hu, L. Reichel, Krylov-subspace methods for the Sylvester equation, 
Linear Algebra Appl. 172 (1992) 283-313. 

[22] I. M. Jaimoukha, E. M. Kasenally, Krylov subspace methods for solving 
large Lyapunov equations, SIAM J. Numer. Anal. 31 (1) (1994) 227-251. 

[23] M. Hochbruck, G. Starke, Preconditioned Krylov subspace methods for 
Lyapunov matrix equations, SIAM J. Matrix Anal. Appl. 16 (1) (1995) 
156-171. 

[24] V. Simoncini, The Lyapunov matrix equation. Matrix analysis from a com¬ 
putational perspective, Tech. Rep. arXiv: 1501.07564, arXiv.org (2015). 
URL http://arxiv.org/abs/1501.07564 

[25] Y. Saad, M. H. Schultz, GMRES: A generalized minimal residual algorithm 
for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 7 
(1986) 856-869. 

[26] H. V. der Vorst, BI-CGSTAB: A fast and smoothly converging variant of 
BI-CG for the solution of nonsymmetric linear systems, SIAM J. Sci. and 
Stat. Comput. 13 (2) (1992) 631-644. 


16 



[27] F. De Teran, F. M. Dopico, Consistency and efficient solution of the 
Sylvester equation for ^-congruence, Electron. J. Linear Algebra 22 (2011) 
849-863. 

[28] V. Simoncini, D. B. Szyld, Theory of inexact Krylov subspace methods and 
applications to scientific computing, SIAM J. Sci. Comput. 25 (2) (2003) 
454-477. 

[29] D. Kressner, C. Schroder, D. S. Watkins, Implicit QR algorithms for palin¬ 
dromic and even eigenvalue problems, Numer. Algorithms 51 (2) (2009) 
209-238. 

[30] Y. Saad, Iterative methods for sparse linear systems, SIAM, 1996. 

[31] J. Wu, Theory and Applications of Partial Functional Differential Equa¬ 
tions, Applied Mathematical Sciences. 119. New York, NY: Springer., 1996. 

[32] B. Vandereycken, S. Vandewalle, A Riemannian optimization approach 
for computing low-rank solutions of Lyapunov equations, SIAM J. Matrix 
Anal. Appl. 31 (5) (2010) 2553-2579. 

[33] L. Hogben (Ed.), Handbook of linear algebra, 2nd Edition, Discrete Math¬ 
ematics and its Applications (Boca Raton), CRC Press, Boca Raton, FL, 
2014. 


Appendix A. Technical bounds 

The following results are needed in the proof of Theorem [7] 

Lemma 8. Suppose Z\ and Z 2 satisfy ([5]) with initial condition Zi(0) 


X. For i = 1,2, 


Z 2 (0) = 


\\Zi(t)\\ F < 2exp(2t(||A 0 || + ||A 1 ||))|| X\\ F . 


Proof. We rely on the vectorized form © of the ODE defining Zi(t ); we have 


II Mil F — 


vec Zi(t) 
vec Z 2 {t) T 


< 11 exp (LA) || 


vec A 
vecX T 


< 2exp(t||A||)||A|| F . 


To complete the proof, we have to estimate the norm of the matrix A in ([T]): 
we have 


||-A|| fs HA? <8 I n \\ + ||(8 drill + ||dn ® Af || + \\I n ® Aq || — 

2(|| A 0 1| + || A, ||), 

where we have used the fact that \\M ® A|| = ||M||||iV||. □ 

Lemma 9. Suppose that Aq has no Hamiltonian eigenpairing. Then, there 
exists a constant K depending only on Aq and c such that 

\\L-\Z)\\ F <KeMr\\A 0 \\/2)\\Z\\ F . 
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Proof. Under the stated hypotheses, T is invertible. Let K be the operator 
norm of T _1 , i.e., the smallest constant such that ||T _1 (Z)||f < I\ ||Z||f- Then 

(A.l) || L~\Z)\\ f = \\T-\Z)eMrA 0 /2)\\ F < 

\\T-\Z)\\ f \\ exp(rA 0 /2)|| < K\\Z\\ F exp(r||A 0 ||/2), 

where we have used the mixed matrix norm inequality ||M./V r || i? < ||M|| F ||./V|| 
123 Page 50-5, Fact 10], □ 
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